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Summary 

We present the development of extended diffraction tomography, a new approach to the solution 
of the linear seismic waveform inversion problem. This method has several appealing features, 
such as the use of arbitrary depth-dependent reference models and the decomposition of the full 
2D or 3D inverse problem into a large number of independent ID problems. This decomposition 
makes the method naturally highly parallelizable. Careful implementation yields significant 
robustness with respect to noise. Several synthetic examples are shown which characterize the 
benefits of our method and demonstrate the usefulness of choosing realistic ID reference media. 



1 Introduction 

There is a vast amount of information contained in seismic data about the Earth structure 
through which the waves have traveled. Traditional data processing has focused on the imaging 
of discontinuities — migration — and the reconstruction of the smooth component of the wave 
velocity variations through the tomographic analysis of traveltimes. However, the use of the 
complete seismic record to obtain high resolution velocity models through full-waveform inversion 
has long been limited by both the computational cost and by a lack of theoretical guidance on 
this nonlinear problem. In recent years, the rapid increase in computational power has enabled 
the development and implementation of some advances in this area, e.g. gradient techniques like 
those developed by Pratt (1999a, 1999b) and Sirgue and Pratt (2004). However, these methods 
are slow to converge and hard to implement in three dimensions. 

More recently, advances in the area of inverse scattering theory (Weglein et al, 2003; Schlott- 
mann, 2006) have arisen that may have the potential of solving the nonlinear full-waveform 
inversion problem more quickly and efficiently. However, these methods rely on having an efficient 
and reliable method of solving the linear waveform inversion problem, i.e., the problem that 
results from assuming that the wave propagation can be approximated by first-order scattering, 
yielding a linear relationship between the data and the unknown subsurface structure. 

The broad class of techniques that address this linear problem fall under the rubric of "diffrac- 
tion tomography". Although the literature on this subject is vast, our work is chiefly influenced 
by Wu and Toksoz (1987). In that work the authors extend earlier work of Devaney (1984) and 
Devaney and Oristaglio (1984) to linearly invert surface reflection profile (SRP) data assuming 
point sources. Wu and Toksoz essentially treat the linear relationship that maps subsurface 
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structure into data as an integral transform and proceed to develop a partial inverse transform. 
In this paper, we will refer to their result as "classical diffraction tomography" . 

The inverse transform of classical diffraction tomography is partial in the sense that not all of 
the subsurface structure can be recovered under realistic conditions. Specifically, in the absence 
of very low-frequency data — the norm in SRP data — many long wavelength components of the 
structure are irretrievable. The origin of this problem is closely tied to the assumption that 
the true subsurface is close to a homogeneous medium and only differs from such by a small 
enough degree that the linear first-order scattering approximation is valid. This assumption of 
a homogeneous background or reference medium in addition to the lack of resolution of long 
wavelengths greatly limits the applicability of this method. 

However, classical diffraction tomography does have the advantage of being an efficient and 
stable method. It is desirable to preserve those qualities while seeking a new theory that frees 
us from the assumption of a homogeneous reference and enables significantly better resolvability 
of subsurface structure. 

In this paper, we present the development of such a result. Our method allows for any ID 
(i.e. depth-dependent) reference medium while also decoupling the full 2D or 3D inverse problem 
into a large number of independent ID inverse problems, resulting in an efficient algorithm that 
is naturally highly parallelizable. With careful implementation, the algorithm that is obtained 
from our theory is stable with respect to noise. By assuming even a marginally realistic ID 
reference medium, the ability to resolve long wavelength components of the subsurface without 
unrealistically low frequencies is greatly improved. 



For the moment, we limit ourselves to a seismic surface reflection profile in only two dimensions. 
We begin by assuming that we measure the scattered component of the seismic wavefield over a 
finite region of the surface. We also assume that the data have been organized in the midpoint- 
offset domain. Specifically, if the data were recorded at horizontal shot /receiver locations xs,xr 
(all located at depth z = 0), we order the data in terms of the coordinates m, h: 
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\h\ < b. 



(2) 
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In the case of most seismic surveys, such as a standard marine survey, the range of the midpoint 
m is controlled by the length of the survey line, whereas the range of h is controlled by the 
aperture of the survey, e.g. the length of a streamer cable. 

Since our objective is to develop an accurate, efficient, and stable solution of the linear 
inverse problem, we take the scattered component of our wavefield to be due to first-order Born 
scattering. In the frequency domain, assuming constant density and acoustic propagation, we 
have 

/oo poo 
dx dz G(xr, zr = 0, u\ x, z)V(x, z)G(x, z, u; xs, zs = 0) (3) 
-oo J 

where the scattering potential V is defined to be 

V(x,z)= 1 27~v ( 4 ) 

c 2 {x,z) c£{z) 

The wave velocities c and Cq are those of the actual and reference media, respectively, and the 
reference medium is assumed to vary only with depth. The function G(x, u;;x ) is the Green's 
function of the reference medium and is the solution to 



cu 2 



G = V 2 G + S (x - x ) . (5) 



4{z) 

Because our reference medium is laterally invariant, G has the important and useful property 

G(x, z,u;x ,z ) = G(\x - x \,z,uj;0,z ). (6) 

Note that the inclusion of a flat free surface will not destroy this property. Using this invariance, 
we rewrite eq. ([3]) in the m, h system (with the source/receiver depths suppressed): 

/oo roo 
dx / dz G(m + h — x, z,u)V(x, z)G(m — h — x, z,cu). (7) 
-oo Jo 

This is the form with which we will work. 

If the wavefield is sampled densely enough relative to the frequency range we consider, we 
can treat m and h as continuous variables. Doing so, we take the Fourier transform of the data 
over the given ranges of m and h: 

D(k m , k h ) = -!- r dm [ b dh e- ikmm e- ikhh D{m, h). (8) 

Z7T J -a J-b 

(Note that we use a symmetric normalization of 1/v2tt in our forward and inverse Fourier 
transforms.) We also take 



G(x - x , z, u; 0, 0) = -= / dk e ik{x - xo) G(k, z, u). (9) 

v27r J-oo 
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Using both of the above in eq. ([3]) and performing the integrals over m and h, we get 
u 2 ab r°° , r°° 



D(k m , k h ) 



oo roo 



dx dz dk dk' sine [(k m — k — k')a] sine [(kh — k + k')b] 

IT J —oo JO J —oo J —oo 



xe- ix(k+k,) V(x, z)G(k, z, uj)G(k', z, uj) (10) 



Making the change-of- variables 



ki. 



k + k', 



(ii) 



(which has a Jacobian determinant of -|) and performing the integral over x, we get 



uj 2 ab 



oo roo 



D(k m ,k h ) = r ^ 1 v27r | dz j dk' m j dk' h sine [(k m - k' m )a] sine [(k h - k' h )b] V(k' m , z) 



xG 



G 



■^{k' m -k' h ),z,uj 



(12) 



For our purposes, the first sine function in eq. ffl2|) is very inconvenient because it couples 
together all the horizontal wavenumbers of the unknown potential V, breaking what would 
otherwise be a complete decoupling of the inverse problem for each value of k m . Fortunately, 
seismic surveys tend to be done over fairly large horizontal areas, which means that a can be 
expected to be large (unlike the aperture). As a tends to infinity, the first sine function will 
tend towards a 5-function. However, in many cases the survey length will still not be large 
enough to allow the 5-function approximation, but we can make another modification to rectify 
the situation. Specifically, by adding a small imaginary part to the frequency, 

uj — > uj + ie, (13) 

we can force the wavefield to decay rapidly enough that we can legitimately make the substitution 

asinc [(k m - k' m )a] ir5 (k m - k' m ) , (14) 

in effect treating our survey length as infinite. Of course, the addition of the imaginary part 
to the frequency is equivalent to tapering our data in the time domain with a slowly decaying 
exponential. In practice, the imaginary part will not need to be much more than a few percent 
of the real part of the frequency. 

What we have gained for our efforts so far is that we can now write the relationship between 
our data and the scattering potential as 



where 

M(k h ,z; k ri 



D(k m ,k h )= dz M(k h ,z;k m )V(k m ,z) (15) 
J o 



uj 2 b f°° ~ rl i ~ rl 

—7= j dk' h sine \{k h - k' h )b] G -(k m + k' h ), z,uj G -(k m - k h ), z,uj . (16) 

V Z7T J— oo IZ J IZ 
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Thus, what we have achieved is the transformation of the 2D inverse problem into a set of 
independent ID inverse problems, one for each horizontal wavenumber k m . 

Before we can use our results, we must consider an additional issue. The symmetry of the 
reference Green's function G with respect to the horizontal offset x — Xq, i.e., its dependence on 
the absolute value \x — xq\, yields the same symmetry in the horizontal wavenumber domain: 

G(k,z,u) = G(-k,z,u). (17) 

This property, after a little thought, implies that M in eq. ffl6|) is symmetric in k m : 

M(k h , z- k m ) = M(k h , z\ -k m ). (18) 

Because of this and because the data have no particular symmetry with respect to k m , if we 
were simply to use eq. ( lT5i) . the appropriate symmetry of the function V, which is the Fourier 
transform of a real-valued function, would not be preserved. In other words, we must reorganize 
the problem to guarantee that 

V(k m ,z) = V*(-k m ,z). (19) 



To do this, we take eq. (I15jl and itself evaluated at — k m and either add or subtract them to 
get two equations: 

roo _ i _, 

/ dz M(k h ,z;k m )Re{V(k m ,z)} = -{D{k m ,k h ) + D{-k m ,k h )) 
Jo 2 

roo _ I „ _ 

/ dz M{k h ,z;k m )Im{V{k m ,z)} = -(D(k m , k h ) - D(-k m , k h )). 20 

This set of equations can be safely used to solve for V for all nonnegative k m , using eq. (fl9~|) to 
retrieve the results for negative wavenumbers. 



2.2 A practical 3D geometry 

Ideally, the extension of our results in 2D to three dimensions would assume that for each shot, 
the wavefield would be recorded over a 2D grid on the surface. However, it is not usually practical 
to cover a large 2D surface area with a dense set of receivers, so we consider instead a 3D survey 
comprised of many 2D surveys. Specifically, we will assume that we have a large number of data 
sets like those we considered in the 2D case, where the sources and receivers are aligned along 
the a>direction, but each set will be for a different value of the y-coordinate with the limitation 
\y\ < c. Thus, the new version of eq. ([3]) is 

/oo roc roo 

dx / dy' / dz' G(x R , y, 0, w\ x\ y', z')V(x', y', z')G(x', y' , z', cj; x s , 0). 
-oo J —oo JO 

(21) 

We continue to use a ID reference model, Cq(z), so G now has the property 

G(x,y,z,u;x Q ,y Q ,Zo) = G(\x - x \, \y-y \,z,w;0,0,zo). (22) 
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We proceed as in the 2D case and organize each survey line into the m, h domain. We again 
take the Fourier transform, 



D[k m , kh, ky 



(27T) 



dm dh dy e- ikmm e- ikhh e- ik « y D(m,h,y), 

J—b J —c 



and expand G, 



1 roo roc _ 

G(x,y,z,cu) = — / dk' x / dk' y e ik * x e ik vVG(k' x ,k' y ,z,cu). 

Z7T J— oo J —oo 



(23) 



(24) 



We also repeat our trick of adding an imaginary part to the frequency so that we can essentially 
take both the survey lengths in both horizontal directions to be infinite. The result, after 
implementing a double set of change-of- variables of the form in eq. ( iTTjl . is 



D[k m , kh, ky) 



dk' 



xG 
xG 



dk'l sine [(k h - k' h )b] V(k m , k y , z') 



(-(^m + k' m ), -{k y + k' h ), z') 



(25) 



As in the 2D case, we have achieved a decoupling of the inverse problem. Here we have broken 
down the full 3D problem into a large set of independent ID problems, one for each pair of 
k m , k y values. Also as with the 2D result, care must be taken to ensure that we account for the 
symmetries in the problem, similar to the procedure used to obtain eqs. ([20]) . 

As a final thought on this 3D treatment, we suggest that improved results would almost 
certainly be obtained by performing an additional set of 2D surveys oriented perpendicular to 
the first. These surveys, aligned along the y-direction, would provide a second, independent set 
of equations for V, most likely improving the condition of the inverse problem. 



3 Numerical Examples 

3.1 Acquisition parameters and methodology 

In this section, we present several simple numerical examples designed to explore the accuracy 
and stability of extended diffraction tomography. In all cases we limit ourselves to the 2D case. 

In each case, the synthetic data is generated via first-order Born scattering (eq. [7j) directly 
in the m, h domain for frequencies of 1-10 Hz at a spacing of 1 Hz. For all our examples, we 
have used the values a = 20 km and 6 = 3 km, indicating a total survey length of 40 km and an 
aperture of 6 km, both of which are reasonable values for modern marine acquisition geometries. 
We have assumed that m and h are regularly spaced at 50 m intervals. All sources and receivers 
are taken to lie at depth z = km. 
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As mentioned in the last section, we add a small imaginary component to the frequency. In all 
cases, that component was .15 s" 1 , which implies an exponential taper in the time domain which 
drops by a factor of e -1 in a characteristic time of 6.67 s. Note that this value was determined 
by trial-and-error to be approximately the smallest effective one. 

Although extended diffraction tomography can incorporate a flat free surface in its reference 
model, we did not do so in these examples but will include that additional bit of realism in future 
work. We used two simple reference models to demonstrate the advantage of being able to have 
a non-constant reference. We begin with a homogeneous medium of constant velocity Co = 1500 
m/s, i.e., water. We follow that with a two-layer medium which has water everywhere above 
depth z = 1 km and C\ = 4000 m/s below. 



To obtain our inversions, we discretized eqs. ([20]) . using the values of k m and kh that naturally 
result from discretizing m and h. The depth coordinate was chosen to range from z = 800 — 
4000 m with a grid spacing of 50 m. Note that the omission of z values less than 800 m 
represents a priori information about V. Although the dependence on frequency was suppressed 
in eqs. (120"]) . we did use multiple frequencies simultaneously in each inversion, a procedure which, 
not surprisingly, augments both stability and accuracy. However, the exact subset used was not 
always the same, as will be discussed later. 

In most of our examples, we also demonstrate the usefulness of filtering in the vertical 
wavenumber (or k z ) domain. The idea is to limit the inversion so that the 2D Fourier transform 
of V, V(k x , kg), is non-zero only within some disk around k = 0. This is done by rewriting the 
integrals over z in eqs. (|20]) as integrals over k z and substituting M and V with their Fourier 
transforms over z. This is another imposition of an a priori bias on our solution, one which helps 
to augment the low-wavenumber part of the solution, which tends to be the most difficult to 
reconstruct. In our examples, the disk was limited to 

|k| < 2.2^ (26) 

where 0J max was the largest angular frequency used in the inversion. Of course, this choice 
was motivated by classical diffraction tomography, where the largest value of |k| obtainable is 
2coWr/co, but in our case we expect to recover some of the "energy" in V(k x ,k z ) outside this 
envelope, so the limit was chosen to be 20% larger. 

Finally, we must mention that the numerical method we used to solve the linear eqs. (120]) was 
the well-known technique of singular value decomposition, or SVD (see e.g. Press et ai, 1992). 
SVD uses the fact that any M x A matrix A such that M > N may be decomposed into the 
product 

A = U S V T (27) 

where U is an M x N column-orthogonal matrix, V is an N x N orthogonal matrix, and S is 
an A x A diagonal matrix whose elements, known as the singular values, are all nonnegative. 
With this decomposition, a solution to the problem 



Ax = d 



(28) 
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may be obtained in a least squares sense: 



x = V S 1 U T d 



(29) 



Since the singular values may be zero or very small, in practice one actually replaces S 1 with 
another matrix S _1 such that if Sj is the j-th singular value, then 



for some tolerance a. 

The selection of the tolerance is quite critical in our inversion problems. Generally, the matrix 
equations that we must solve are very poorly conditioned, with many singular values being very 
close to zero. The main difficulty that arises from this poor conditioning is an extreme sensitivity 
to noise. This occurs because for very small Sj, the j-th diagonal element of S _1 is extremely 
large, amplifying any errors in our data, such as noise. In our computations, a tolerance of a — 1 
was used uniformly. As we will see, this will provide robustness with respect to noise. However, 
the price for this robustness is some degradation in the quality of our inversion, as some features 
of V which we would like to recover only weakly affect the data through some of the small sj 
that we drop from consideration. 



3.2 Examples 

In all of the following figures, we will see images of both the true models used and the inversion 
results as well cross-sections through them. In the model and inversion images, only the middle 
4 km of the full 40 km re-coordinate range is shown in order to allow the inversion results to be 
more easily visible. Although it is not indicated in the figures, as stated above, the sources and 
receivers are assumed to be located at z — 0, the top of the image boxes. Each cross-section 
shown, vertical and horizontal, is taken through the center of the anomaly we are trying to 
reconstruct. Note that in each such slice, the true model will always be plotted in red while the 
inversion results will be in blue. 

In all cases, the input model is a square-shaped anomaly of side length 1 km centered hor- 
izontally at x = but with different depth for different examples. The value of the scattering 
potential V, which is a dimensionless quantity, is always 1, corresponding to a low- velocity region. 
However, since this is a purely linear inverse problem, the actual magnitude and sign of V are 
irrelevant. This would most certainly not be the case for the nonlinear full- waveform problem. 

We begin by looking at some results for a homogeneous reference medium. In Figures 1 and 
2, we see the results for our square centered at a depth of z = 1500 m. In this test, frequencies 
of 1-5 Hz were used. Notice that there is marked improvement in the results when the range of 
k z values is restricted. 




(30) 
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Figure 1: Results for a homogeneous background using frequencies 1,2,. ..,5 Hz. (a) True model, 
(b) Result for unrestricted k z . (c) Result for restricted k z . 



(a) (b) 




1 2 3 4 -2 -1 1 

z (km) x (km) 



Figure 2: Cross-sections through the centers of the squares in Fig. 1. In each case, the true 
model is in red, and the reconstruction is in blue. (a,b) Vertical and horizontal slices through 
the unrestricted result. (c,d) Vertical and horizontal slices through the restricted results. 
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In Figures 3 and 4 we have the results for the same homogeneous background but with 
the square now centered 500 m deeper at z — 2000 m. Clearly, the inversion suffers serious 
degradation from even this slight increase in target depth. We are already reaching the limit of 
our method for a constant velocity reference, even with k z restriction. These images could be 
improved by lowering the tolerance a, but then sensitivity to noise would grow rapidly. 

Next, in Figures 5 and 6, we explore the effect of having a two-layer reference model on 
reconstructing the deeper anomaly. Note that Figure 5a shows the scattering potential V, while 
Figure 5b shows the complete velocity model. With nothing else being changed from the previous 
example, we see dramatic improvements in the image quality. We attribute this improvement 
to two effects. One is the "stretching" of the wavelengths of the incident waves as they enter 
the higher-velocity lower half-space. Thus, for the same frequency you gain sensitivity to lower- 
wavelength components of V. (This actually results in slightly worse reconstruction of the 
higher wavelengths of V '.) The other effect which adds to our image quality is the refraction of 
the incident waves, which brings them closer to horizontal, again improving our sensitivity to 
wavelength components that are hard to recover in a homogeneous reference model. Note also 
that the restriction of k z adds little to the improvement of the results here. 

In Figures 7 and 8, we see the advantage of adding higher frequencies into the inversion. Here 
we have used the full range of frequencies, 1-10 Hz, that were prepared for these tests. Clearly, 
the anomaly is well-resolved at all wavelengths, and, again, k z restriction adds little other than 
the dampening of small artifacts. 

In real seismic exploration data, we are usually limited to frequencies larger than 5 Hz. In 
Figures 9 and 10, we see the results of extended diffraction tomography when we use only 5-10 
Hz frequencies and the same two-layer model as before. We do see degradation of the results, but 
k z restriction helps to subdue them partially. We hypothesize that using an even more realistic 
ID reference model, one with a gradient in the lower half-space, would significantly improve 
the results because of the presence of turning rays, which would sample yet more effectively the 
components of the spectrum of V that are presently most difficult to obtain. Future work will 
explore this hypothesis. 

Our last example demonstrates the robustness of our method. Here we test the effect of 
adding 5% noise. By that, we mean that for each frequency, we added to D(m,h;uj) complex- 
valued numbers that had a constant absolute value of .05 times the maximum absolute value of 
D at that frequency but had a phase that was uniformly randomly distributed between and 
2vr. 

In Figure 11, we see the effect of adding the noise to our data for the two-layer reference 
model. In this case, we used the full 1-10 Hz frequency range. Although there is now noise in 
the recovered image, our use of a large tolerance in the SVD inversion has preserved the vast 
majority of the structure. 
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Figure 3: Results for the same situation as in Fig. 1 except that the square is 500 m lower, (a) 
True model. (b,c) Unrestricted and restricted results. 
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Figure 4: Cross-sections through the results of Fig. 3. (a,b) Vertical and horizontal slices through 
the unrestricted result. (c,d) Vertical and horizontal slices through the restricted results. 
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Figure 5: Results for the two-layer background model. Again frequencies 1,2,.. .,5 Hz are used, 
(a) True model of V. (b) True velocity model (V plus background). (c,d) Results for unrestricted 
and restricted cases. 



(a) (b) 
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Figure 6: Cross-sections through the results of Fig. 5. (a,b) Vertical and horizontal slices through 
the unrestricted result. (c,d) Vertical and horizontal slices through the restricted results. 
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Figure 7: Results for the two-layer background model but with frequencies 1,2,. ..,10 Hz used, 
(a) True model of V. (b,c) Results for unrestricted and restricted cases. 




Figure 8: Cross-sections through the results of Fig. 7. (a,b) Vertical and horizontal slices through 
the unrestricted result. (c,d) Vertical and horizontal slices through the restricted results. 
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Figure 10: Cross-sections through the results of Fig. 9. (a,b) Vertical and horizontal slices 
through the unrestricted result. (c,d) Vertical and horizontal slices through the restricted results. 
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Figure 11: Results for the two-layer background model for 1,2,. ..,10 Hz and 5% noise added, (a) 
True model of V. (b) Reconstruction result. (c,d) Vertical and horizontal slices through the 
result. 



4 Discussion 



We have presented the development of extended diffraction tomography, a new variant of diffrac- 
tion tomography that can very efficiently and often very accurately solve the linear seismic 
inversion problem for surface reflection data with finite aperture. The algorithm, when properly 
implemented, exhibits excellent robustness with respect to random noise. In addition, we have 
extended the theory to a realistic 3D acquisition geometry. 

The main strength of extended diffraction geometry originates from the use of laterally invari- 
ant (i.e., ID) reference models. This allows the decomposition of the full 2D inverse problem into 
many independent ID problems, making the method naturally highly parallelizable. In contrast 
to classical diffraction tomography, the use of ID reference models also allows for greater realism 
and accuracy in the inversion. We have already demonstrated that a simple two-layer reference 
model greatly enhances the results. Further work will be done to confirm that even more realistic 
models containing gradients will yield yet more substantial improvements. 

The reduction of the inverse problem to the solution of relatively small matrix equations opens 
the possibility of applying the formidable tools of numerical linear algebra. We have already 
demonstrated that the simple restriction of the allowed values of the vertical wavenumber k z 
yields significant benefits. Future efforts may include the incorporation of more sophisticated 
methods such as Tikhonov and total variation regularizations (see e.g. Aster et al, 2005). 
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Further testing of the method must also be done to explore how it handles data that is not 
generated via first-order Born scattering, i.e., we must see how it handles multiples, in order to 
evaluate its usefulness in nonlinear inversion. Although this method may have some value in and 
of itself, its ultimate purpose is to be used in the full-waveform nonlinear inversion of realistic 
seismic data. 
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